A blue sky catastrophe in double-diffusive convection 
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A global bifurcation of the blue sky catastrophe type has been found in a small Prandtl number 
binary mixture contained in a laterally heated cavity. The system has been studied numerically 
applying the tools of bifurcation theory. The catastrophe corresponds to the destruction of an orbit 
which, for a large range of Rayleigh numbers, is the only stable solution. This orbit is born in 
a global saddle-loop bifurcation and becomes chaotic in a period doubling cascade just before its 
disappearance at the blue sky catastrophe. 

PACS numbers: 47.27.Te, 47.20.Ky, 44.25.+f 



Bifurcation theory has long been a very helpful tool in 
the analysis of complex dynamics of nonlinear systems 
0, Whereas different devised scenarios have been 
found in theoretical models with a few variables, there is 
a growing interest both in relating real systems with that 
kind of models (e.g. projecting their dynamics to some 
relevant degrees of freedom |3() and in directly analyz- 
ing the behavior of these systems in terms of dynamical 
systems theory (by studying them either experimentally 
or by realistic models). In this context a great deal of 
work has been devoted to convection in fluids. Qualita- 
tive changes in the dynamics of fluxes maintained out of 
equilibrium by imposed thermal gradients have provided 
examples of most of the known bifurcations, and have be- 
come a main subject in the area of nonlinear dynamics. 

In this letter we will show the occurrence of a blue sky 
catastrophe [BSC] in double diffusive convection. The 
BSC is a codimension-1 bifurcation that consists in the 
destruction of a stable periodic orbit as its length and 
period tend to infinity, while the cycle remains bounded 
and located at a finite distance from all the equilibrium 
solutions 0, 0| . This destruction is caused by the col- 
lision with a non-hyperbolic cycle that appears at the 
bifurcation point. While approaching the bifurcation the 
orbit increasingly coils in the zone where the new cycle 
will appear, which originates the divergence in both pe- 
riod and length. In that point the original cycle becomes 
an orbit homoclinic to the new cycle. This type of bi- 
furcation is relatively exotic, but can easily be found in 
slow-fast (i.e. singularly perturbed) systems with at least 
two fast variables 

We are interested in double-diffusive fluxes that occur 
when convection is driven by simultaneous thermal and 
concentration gradients in a binary mixture . Double- 
diffusive convection in cavities with imposed vertical gra- 
dients exhibits very rich dynamics, and has been used 
as a system to study pattern formation and transi- 
tion to chaos || . The case of horizontal gradients, which 
arises naturally in applications such as crystal growth 
or oceanography has received less attention. In this 



work we numerically study this latter configuration for 
a small Prandtl number binary mixture. We consider 
the case when thermal and solutal buoyancy forces ex- 
actly compensate each other, which allows the existence 
of a quiescent (conductive) state 0, 0, 0, 0. We 
have found that in this system there exists a large range 
of Rayleigh numbers in which the only stable solution 
is an orbit that features a low-frequency spiking behav- 
ior. This orbit appears associated to a global bifurcation 
and loses stability when a period doubling cascade takes 
place originating a chaotic attractor. However, the most 
remarkable feature of this chaotic attractor is its sudden 
disappearance in a BSC of the chaotic type. As far as we 
know this is the first example of such bifurcation in an 
extended system. 

We have considered a binary mixture in a 2-D rect- 
angular cavity of aspect ratio T = d/h = 2, where d is 
the length and h is the height of the cavity. A difference 
of temperature AT is maintained between both vertical 
boundaries. Dimcnsionless equations in Boussinesq ap- 
proximation explicitly read 
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where u = (v x ,v z ) is the velocity field in (x,z) coordi- 
nates, P is the pressure over the density, 9 denotes the de- 
parture of the temperature from a linear horizontal pro- 
file. C is the scaled deviation of the concentration of the 
heavier component relative to the linear horizontal profile 
which equilibrates that of the temperature in the expres- 
sion of the mass flux. Lengths and times are scaled with h 
and t K — h 2 /k, respectively, being k the thermal diffusiv- 
ity The dimensionless parameters are the Prandtl num- 
ber <7 = v/k, the Rayleigh number Ra = agh 3 AT/vK 
and the Lewis number r = D/n, where v denotes the 
kinematic viscosity, g the gravity level, a the thermal ex- 
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FIG. 1: Nusselt number of steady solutions versus Ra, and its 
corresponding bifurcations. Continuous lines: stable states. 
Dashed lines: unstable states. 



pansion coefficient, and D is the mass diffusivity. The 
separation ratio S = Cq(1 — Cq)^St will be taken 
S = — 1. Here, St is the Soret coefficient, C is the 
concentration of the heavier component in the homoge- 
neous mixture, and is the mass expansion coefficient 
(j3 > for the heavier component). 

The boundaries are taken to be no-slip and with no 
mass flux. Lateral walls are maintained at constant tem- 
peratures and at the horizontal plates a linear profile of 
temperature between the two prescribed temperatures is 
imposed. Thus, boundary conditions are written as 



n-V(C-0)=O, at dn. 



(5) 



Notice that these boundary conditions prevent one to 
absorb the Soret terms into the equations like in Refs. 
[HI El El E3 On the other hand this system is 
Z2-equi variant. Eqs. 1114(1 . together with boundary 
conditions JSJ), are invariant under a transformation 7r, 
a central symmetry around the point (r/2,1/2), i.e. 
ir : (v x ,v z ,6,C) -> (-v x , -v z , -6, -C), (x, z) -» (T - 
x, 1 — z). Hence any solution of these equations either is 
7r-invariant (for now on we will call it symmetric) or its 
image under 7r is also a solution (constituting a pair of 
asymmetric solutions). This has important consequences 
on the nature of its possible bifurcations pj. 

We have obtained time-dependent solutions of equa- 
tions 1114(1 and boundary conditions J5J by using a second 
order time-splitting algorithm, proposed in Ref. 1 1 ■ ll ] , and 
a pseudo-spectral Chebyshev method for the space dis- 
cretization. Furthermore we have calculated (both sta- 
ble and unstable) steady solutions and analyzed their 
stability by adapting a pseudospectral first-order time- 
stepping formulation, as described in Ref. [lil Hfil Il7j. 
The values of the parameters have been a = 0.00715 and 
t = 0.03, close to that characteristic of molten doped 
germanium [lil fl9j| . Spatial discretization has typically 
been between 60 x 30 and 90 x 60 mesh grid points. 

The scenario provided by the analysis of the steady 





FIG. 2: Concentration levels of the steady solutions of the 
symmetric branch (Ra = 1252) and the non symmetric branch 
(Ra ■- 



solutions is shown in the bifurcations diagram of Fig. ^ 
In this figure the Nusselt number (Nu), defined as the 
quotient of heat flux through the hot wall with that of 
the corresponding conductive solution, is represented for 
the steady states as a function of the Rayleigh number 
(Ra). For the sake of clarity only one asymmetric so- 
lution of each pair has been shown. For small Ra the 
conductive solution (allowed here by the choice S = — 1) 
is stable, but loses stability, maintaining the symmetry, 
through a transcritical bifurcation at Ra — 541.9. The 
supercritical branch of the bifurcating solution is only 
stable up to a pitchfork bifurcation at Ra = 542.4, fol- 
lowing a scenario similar to that described in Ref. |l2| . 
The interesting behavior in this system originates from 
the subcritical branch. This branch gains stability via a 
saddle node bifurcation at Ra — 99 (SNi), and loses it 
again at Ra = 245 in a Pitchfork bifurcation (P) where 
a couple of stable asymmetric branches appear. In Fig. 
121 we represent the concentration for a symmetric (left) 
and an asymmetric (right) steady states. We can see 
that concentration is roughly homogeneous inside rolls, 
displacing concentration gradients to the lateral bound- 
aries. 

The asymmetrical steady state is stable until Ra = 
1209, where it loses stability at a saddle-node bifurcation 
(S7V2). The full branch of asymmetrical steady states 
is depicted in Fig. ^ where we can see that it changes 
again the direction at a turning point at Ra = 865.6, 
but without gaining stability. Increasing the Rayleigh 
number Hopf bifurcations of the symmetric and asym- 
metric branches take place at Hi (Ra = 2137) and H2 
(Ra = 2218) respectively. The branch of symmetric peri- 
odic orbits emanating from H± will play an essential role 
in the subsequent evolution of the system. 

In the range from Ra = 1209 until Ra — 2253 we 
have found no stable solution connected with the above 
branches by local bifurcations. Integrating the evolu- 
tion equations we have obtained a branch of asymmetric 
periodic solutions that dominates the dynamics of the 
system in this range of parameters. In Fig. [3] we rep- 
resent time series and phase space plots of the orbits of 
this branch for two different values of the Rayleigh num- 
ber. The oscillations first appear in the form of spikes of 
very large period (see Fig. [3| a), according to the prox- 
imity to a global saddle-loop bifurcation that occurs at 
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FIG. 3: Velocity components of a representative point. Top: 
Asymmetric orbit at Ra = 1183.68. a) time series, with the 
value of the saddle stationary solution marked, b) orbit in 
the phase space. Bottom: Attractor at Ra = 2255. c) time 
series, d) attractor in the phase space with the stable sym- 
metric orbit. The unstable stationary symmetric solution is 
also shown. 
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FIG. 4: Left: logarithmic fit of the periods for the SL con- 
nection. Right: Square-root fit of the periods for the BSC. 



frequency of the unstable symmetric orbit that appears 
in Hi . In fact we have been able to calculate this unsta- 
ble branch by temporal evolution forcing the symmetry 
of the system, and its frequency coincides with that of 
the windings of the attractor on all the branch. If we in- 
crease further Ra, the asymmetric orbit follows a period- 
doubling cascade, becoming chaotic. This is revealed in 
the phase of the winding of the trajectory, as can be seen 
in Fig. OD where a detail of the orbit during the cascade 
is shown. This cascade seems to move to slightly higher 
Ra values as spatial resolution is increased, but we have 
not been able to obtain the precise values due to the 
extremely large duration of the orbits in this regime. In 
Fig. El(c,d) the attractor thus generated is represented at 
Ra = 2255. For this value of Ra the symmetric orbit has 
already become stable at a Pitchfork bifurcation (P so , at 
Ra = 2253), and both coexist. Very shortly afterwards, 
the whole attractor disappears at Ra — 2257.5. 

This destruction of the attractor exhibits characteris- 
tics that permit to identify it as the chaotic counterpart 
of the scenario for BSC bifurcation described in Refs. 
P, |2(J ■ Indeed in all the process the attractor remains 
bounded and at a finite distance of any steady solution, 
as required 0. The average length and time between 
spikes (which are reproducible with variations smaller 
than 1 over 1000) diverges as the windings start to accu- 
mulate, which occurs at a specific location in the attrac- 
tor. That indicates that the solution is colliding there 
with a new cycle that appears at the bifurcation point, 
and to which it becomes homoclinic. Furthermore this 
divergence, shown in Fig. 0] (right), is very well fitted by 
a square-root law: 



Ra = 1183.67 (SL) where the orbit connects with the 
unstable branch of SN2 (see Fig. QJ. The character of 
this global bifurcation can be inferred from the logarith- 
mic divergence of the period when the Rayleigh number 
decreases toward (SL). We have fitted that period to 



T 



A 



log(Ra- Ra SL ) + A. 



(6) 



We can see the fit in Fig. 0] (left). The resulting value 
A/it = 0.079 results to be quite close to the unstable 
eigenvalue A = 0.074 of the saddle stationary point, as 
obtained by the stability calculation. Near that global 
bifurcation the time evolution of the velocity of a repre- 
sentative point is shown in Fig. [3] (a). The value for the 
saddle asymmetric state is also represented. We can see 
how the solution spends a long time near it. The spike 
corresponds to a rapid and large excursion by the phase 
space, as seen in Fig. 01(b), during which the roll alter- 
nately switches between a confined and a more centered 
positions (analogous to the patterns shown in Fig. 0) . 

Increasing Ra, at Ra = 2137 the orbit starts to curl, 
showing ripples in the time dependence, reflecting the 



\J Ra c — Ra 



+ B. 



(7) 



This law of divergence particularly corresponds to that 
scenario, since it demonstrates that the new cycle to 
which the attractor is connecting is the saddle-node of 
two orbits (SNOi). In principle there are several possi- 
bilities for the topology of the attractor . In our case 
the successive windings are braided by the tip of the at- 
tractor into an almost one dimensional tube or filament 
(Fig. |5J). This filament reintroduces the orbit into the 
vicinity of the saddle- node orbit, and it starts winding 
again accumulating curls near it. Therefore, in the limit, 
the attractor has the topology of a French Horn. This 
feature is also shared with Ref. |20| . 

After the BSC, one would expect the system to reach 
the stable member of the pair of asymmetric solutions 
born at SNOi. On the contrary, simulations show that 
the system evolves through a extremely long transient, 
during which the trajectory accumulates curls near the 
saddle-node before being rejected to the symmetric orbit 
that became stable at Pso ■ That could mean either that 
the stability range of the asymmetric orbit is very small 
(which would require a much finer exploration in Ra to 
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FIG. 5: Period-doubling cascade (zoom of the tip of the at- 
tractor). a) period 1 (Ra = 2220). b) period 2 (Ra = 2232). 
c) period 4 (Ra = 2235). d) chaotic solution (Ra = 2240). 
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FIG. 6: Diagram of the conjectured unstable asymmetric or- 
bit (thick line) and its connections to other branches 

find it, a formidable task in this slow regime), or that 
its basin of attraction is very reduced (and the nearby 
symmetric orbit attracted all the calculated orbits). 

We propose that the SNOi is located in the branch 
of unstable asymmetric orbits created at the pitchfork 
bifurcation where the new orbit becomes stable (P so )- 
This hypothetical scenario is shown in Fig. and is 
the simplest one in which the attractor presents at the 
BSC an homoclinic connection to a branch coming from 
known solutions. This conjecture requires the unstable 
asymmetric branch to gain stability in a first saddle node 
bifurcation SW0 2 an d to lose it again at the SNOi, as 
can be seen in Fig. [S] The coincidence of the frequency 
value of the symmetric orbit at Psoi ^so — 7.01, to that 
of the windings of the attractor, uj = 7.01, is consistent 
with such connection. A small distance between SNOi 
and SNO2 would explain the reduced stability domain 
of the asymmetrical solution. Finally, the presence of 
the additional saddle node orbit SNO2 in the proximity 
would slow down the dynamics for Ra slightly above, 
making the transient to the symmetric solution very long, 
as it is actually observed. 

One could devise more complex scenarios for the oc- 
currence of this BSC. For example, the attractor could be 



destroyed on a boundary crisis associated to global con- 
nections of the unstable orbits coming from the period 
doubling bifurcations. 

Finally, it is worth remarking than the BSC displayed 
by this system is robust against small changes in the value 
of the separation ratio S. In particular we have obtained 
a similar BSC in simulations performed with S = —0.99. 
That means that the additional symmetry introduced in 
the system by the special value S — — 1 is not an essential 
ingredient of the phenomena described here. 
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